{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "eight-monkey", "metadata": {}, "outputs": [], "source": [ "# This tutorial will use Monte Carlo methods to calculate the volume of an\n", "# n-dimensional sphere.\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "id": "ideal-taste", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.7581937008533004" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# In all Monte Carlo simulations it is necessary to generate random or\n", "# pseudo-random numbers. The Python command 'np.random.uniform()'\n", "# from the NumPy module generates random numbers uniformly distributed \n", "# between zero and one. \n", "np.random.uniform()" ] }, { "cell_type": "code", "execution_count": 3, "id": "equal-somewhere", "metadata": {}, "outputs": [], "source": [ "# This tutorial will attempt to numerically determine the volume of a sphere\n", "# of arbitrary dimension using the Hit & Miss method. This is equivalent to\n", "# evaluating a d-dimensional integral. In d-dimensions, a hypersphere is\n", "# defined by x1^2+x^2+x3^3+...+xd^2 < R^2.\n", "\n", "# To find the volume of a sphere we'll generate random numbers over some\n", "# \"square\" volume and count how many land inside the sphere. The\n", "# probability of the randomly generated point landing within the sphere is\n", "# simply p = Vsphere / Vsquare. If we use spheres of radius 1, then the\n", "# volume of the sqaure that contains the sphere is 2^d where d is the\n", "# number of dimensions that we're working with. The Monte Carlo simulation\n", "# is used to determine p = Zn/n where n is the number of trials and Zn is\n", "# the number of \"Hits\". Collecting all of these results we have\n", "# p = Zn/n = Vsphere/2^d or Vsphere = 2^d*(Zn/n)." ] }, { "cell_type": "code", "execution_count": 4, "id": "suspected-torture", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The volume of a 1 D sphere from the Monte Carlo simulation is 2.0\n", "The exact volume is 2 and the ratio of the two is 1.0\n" ] } ], "source": [ "# We'll start with the trivial example of a 1-D sphere. In 1-D, x1^2" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Below we import and display an image from Wikipedia that shows the\n", "# surface areas and volumes of \"n-spheres\" up to n = 9\n", "# (https://en.wikipedia.org/wiki/N-sphere).\n", "import matplotlib.image as mpimg\n", "img = mpimg.imread('hypersphere.png')\n", "plt.figure(figsize = (12, 7))\n", "imgplot = plt.imshow(img)" ] }, { "cell_type": "code", "execution_count": 11, "id": "fatal-exclusive", "metadata": {}, "outputs": [], "source": [ "# Now let's do something more interesting... Let's calculate the volume of\n", "# a high-dimensional sphere -- one that we don't know ahead of time the\n", "# correct answer. Let's arbitrarily take d = 12.\n", "d = 12 # The number of dimensions\n", "n = int(1e6) # The number of Monte Carlo trials\n", "Zn = 0\n", "for i in range(n):\n", " x = np.random.uniform(0, 1, d)\n", " Rsq = 0\n", " for j in range(d):\n", " Rsq = Rsq + x[j]**2\n", " R = np.sqrt(Rsq) # Calculate distance from the orgin\n", " if R <= 1: # Check for a hit\n", " Zn += 1 # If there's a hit increment Zn" ] }, { "cell_type": "code", "execution_count": 12, "id": "forbidden-luxembourg", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The volume of a 12 D sphere from the Monte Carlo simulation is 1.339392\n" ] } ], "source": [ "# Here's the result! The volume of this sphere is this prefactor times the\n", "# radius of the sphere to the 12th power.\n", "V12D = 2**d*Zn/n # The Monte Carlo calculation of the sphere volume\n", "print('The volume of a', d, 'D sphere from the Monte Carlo simulation is', V12D)" ] }, { "cell_type": "code", "execution_count": 13, "id": "modern-shark", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2021-03-11 11:32:50.516899\n", "0%\n", "10%\n", "20%\n", "30%\n", "40%\n", "50%\n", "60%\n", "70%\n", "80%\n", "90%\n", "2021-03-11 13:28:22.524509\n" ] }, { "data": { "text/plain": [ "Text(0, 0.5, 'Counts')" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAWwklEQVR4nO3dfbRddX3n8fdHAcGnCpJArMRoRWqrAnqxtfgARlq0jKAjouNDVJy0tdo6jlacznJqZ9YMHV3WWW2VybKa2HFQpGAorWiaAdGKSnhQVMBY5WkZk4i1PlYJfuePs1MONzf3nntz97n35vd+rXXXPvv5e3ZOPmef39n7d1JVSJLacZ+FLkCSNF4GvyQ1xuCXpMYY/JLUGINfkhpzwEIXMIrDDz+8Vq1atdBlSNKScs0113y7qpZNnr4kgn/VqlVs2bJlocuQpCUlya1TTbepR5IaY/BLUmMMfklqTG/Bn+SYJNcP/X0vyeuTHJZkU5Kt3fDQvmqQJO2pt+Cvqpur6riqOg54EvAj4GLgHGBzVR0NbO7GJUljMq6mntXAP1bVrcDpwIZu+gbgjDHVIElifMH/IuD87vERVbUNoBsun2qFJGuTbEmyZefOnWMqU5L2f70Hf5KDgOcCH5nNelW1rqomqmpi2bI97j+QJM3ROM74nw1cW1Xbu/HtSVYAdMMdY6hBktQZR/C/mHuaeQAuAdZ0j9cAG8dQg8bgyJUrSTLnvyNXrlzopyA1odcuG5LcHzgF+K2hyecCFyQ5G7gNOLPPGjQ+22+/HS6/fO7rn3zyPFYjaW96Df6q+hHw0EnT7mRwlY8kaQF4564kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSY3oN/iQPSXJhkpuS3JjkKUkOS7IpydZueGifNUiS7q3vM/7/BVxWVb8IHAvcCJwDbK6qo4HN3bgkaUx6C/4kDwaeDvwlQFX9tKq+C5wObOgW2wCc0VcNkqQ99XnG/yhgJ/D+JNcleW+SBwBHVNU2gG64fKqVk6xNsiXJlp07d/ZYpiS1pc/gPwB4IvCeqjoe+CGzaNapqnVVNVFVE8uWLeurRklqTp/BfwdwR1V9rhu/kMEbwfYkKwC64Y4ea5AkTdJb8FfVt4DbkxzTTVoNfAW4BFjTTVsDbOyrBknSng7oefuvAz6Y5CDg68ArGbzZXJDkbOA24Myea5AkDek1+KvqemBiilmr+9yvJGnvvHNXkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEH9LnxJLcA3wfuBnZV1USSw4APA6uAW4AXVtU/9VmHJOke4zjjP7mqjquqiW78HGBzVR0NbO7GJUljshBNPacDG7rHG4AzFqAGSWpW38FfwCeSXJNkbTftiKraBtANl0+1YpK1SbYk2bJz586ey5SkdvTaxg+cWFXfTLIc2JTkplFXrKp1wDqAiYmJ6qtASWpNr2f8VfXNbrgDuBh4MrA9yQqAbrijzxokSffWW/AneUCSB+1+DPw68CXgEmBNt9gaYGNfNWiJOfBAkszp78iVKxe6emnJ6LOp5wjg4iS79/N/q+qyJFcDFyQ5G7gNOLPHGrSU3HUXXH75nFbdfvLJ81yMtP/qLfir6uvAsVNMvxNY3dd+JUnT885dSWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX/sHb/6SRtZ3Xz3SeHjzlzQyz/glqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9585ca4w1ckjd/qTGe8UtSYwx+SWqMwS9JjZl18Cc5NMkT+ihGktS/kYI/yRVJHpzkMOALwPuTvHPEde+b5Lokl3bjhyXZlGRrNzx07uVLkmZr1DP+n6uq7wHPB95fVU8CnjXiur8P3Dg0fg6wuaqOBjZ345KkMRk1+A9IsgJ4IXDpqBtP8nDgN4H3Dk0+HdjQPd4AnDHq9iRJ+27U4H8b8HHga1V1dZJHAVtHWO9dwB8APxuadkRVbQPohsunWjHJ2iRbkmzZuXPniGVKkmYyavBvq6onVNVrAKrq68C0bfxJTgN2VNU1cymsqtZV1URVTSxbtmwum5AkTWHU4P+zEacNOxF4bpJbgA8Bz0zyf4DtXbMR3XDHiDVIkubBtF02JHkK8GvAsiRvGJr1YOC+061bVW8B3tJt5yTgjVX10iRvB9YA53bDjXMtXpI0ezP11XMQ8MBuuQcNTf8e8II57vNc4IIkZwO3AWfOcTuSpDmYNvir6pPAJ5Osr6pb57qTqroCuKJ7fCeweq7bkiTtm1F757xfknXAquF1quqZfRQlLRldl85zdcRRR/Gt226bx4KkmY0a/B8BzmNwPf7d/ZUjLTH70KUz2K2zFsaowb+rqt7TayWSpLEY9XLOv0nymiQrur52Duv67ZEkLTGjnvGv6YZvGppWwKPmtxxJUt9GCv6qemTfhUiSxmOk4E/y8qmmV9UH5rccSVLfRm3qOWHo8cEMrsO/FjD4JWmJGbWp53XD40l+DvirXiqSJPVqrr+5+yPg6PksRJI0HqO28f8Ng6t4YNA522OBC/oqSpLUn1Hb+N8x9HgXcGtV3dFDPZKkno3U1NN11nYTgx46DwV+2mdRkqT+jBT8SV4IfJ5BF8ovBD6XZK7dMkuSFtCoTT1/CJxQVTsAkiwD/h64sK/CJEn9GPWqnvvsDv3OnbNYV5K0iIx6xn9Zko8D53fjZwF/109JkqQ+zfSbu48GjqiqNyV5PvBUIMBVwAfHUJ8kaZ7N1FzzLuD7AFV1UVW9oar+A4Oz/Xf1W5okqQ8zBf+qqvri5IlVtYXBzzBKkpaYmYL/4GnmHTKfhUiSxmOm4L86yb+fPDHJ2cA1062Y5OAkn0/yhSRfTvK2bvphSTYl2doND517+ZKk2Zrpqp7XAxcneQn3BP0EcBDwvBnW/QnwzKr6QZIDgU8n+RjwfGBzVZ2b5BzgHODNc30CkqTZmTb4q2o78GtJTgYe103+26r6fzNtuKoK+EE3emD3V8DpwEnd9A3AFRj8kjQ2o/bHfzlw+Ww3nuS+DD4pPBr4i6r6XJIjqmpbt91tSZbvZd21wFqAlStXznbXkqS96PXu26q6u6qOAx4OPDnJ42ZYZXjddVU1UVUTy5Yt661GSWrNWLpdqKrvMmjSORXYnmQFQDfcsfc1JUnzrbfgT7IsyUO6x4cAz2LQtfMlwJpusTXAxr5qkCTtadS+euZiBbCha+e/D3BBVV2a5Crggu6S0NsYdPUsSRqT3oK/u+P3+Cmm3wms7mu/kqTp2bWy7uXIlStJMqc/SUtDn009WoK23347XD7rK3cHTj55fouR1AvP+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8O+H7FpZ0nTslnk/ZNfKkqbjGb8kNcbgX6RsrpHUF5t6FimbayT1xTN+SWpMb8Gf5Kgklye5McmXk/x+N/2wJJuSbO2Gh/ZVgyRpT32e8e8C/mNVPRb4VeB3k/wScA6wuaqOBjZ345KkMekt+KtqW1Vd2z3+PnAj8PPA6cCGbrENwBl91SBJ2tNY2viTrAKOBz4HHFFV22Dw5gAs38s6a5NsSbJl586d4yhTkprQe/AneSDw18Drq+p7o65XVeuqaqKqJpYtW9ZfgZLUmF6DP8mBDEL/g1V1UTd5e5IV3fwVwI4+a5Ak3VufV/UE+Evgxqp659CsS4A13eM1wMa+apAk7anPG7hOBF4G3JDk+m7afwLOBS5IcjZwG3BmjzVIkibpLfir6tPA3voPWN3XfiVJ0/PO3R7Z345mdOCBc36NHLly5UJXryXKvnp6ZH87mtFdd835NbLd14jmyDN+SWqMwS8tVTYTaY5s6pGWKpuJNEee8UtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGtNb8Cd5X5IdSb40NO2wJJuSbO2Gh/a1f0nS1Po8418PnDpp2jnA5qo6GtjcjUuSxqi34K+qK4HvTJp8OrChe7wBOKOv/UuSpjbuNv4jqmobQDdcvrcFk6xNsiXJlp07d46tQKkJ/lB70xbtj61X1TpgHcDExEQtcDnS/sUfam/auM/4tydZAdANd4x5/5LUvHEH/yXAmu7xGmDjmPcvSc3r83LO84GrgGOS3JHkbOBc4JQkW4FTunFJ0hj11sZfVS/ey6zVfe1TkjQz79yVpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4p3HkypVz7s8kyUKXL0lTWrR99SwG22+/fc79mQBgnyaSFiHP+CWpMQa/pNmxS+clz6YeSbNjl85Lnmf8ktSY/T749+XKHEnaH+33TT37dGWOH0sl7Yf2+zN+SdK9GfyS1BiDX5IaY/BLUmMMfknjsw83f3kD2PzZ76/qkbSI7MPNX+ANYPPFM35JasyCBH+SU5PcnORrSc5ZiBokLUFLsJ+gfe3evY+6x97Uk+S+wF8ApwB3AFcnuaSqvjLuWiQtMUuwn6B97d69j7oX4oz/ycDXqurrVfVT4EPA6QtQhyQ1KVU13h0mLwBOrapXd+MvA36lql47abm1wNpu9Bjg5rEWum8OB7690EUsQh6XqXlcpuZxmdpsjssjqmrZ5IkLcVXPVL2f7fHuU1XrgHX9lzP/kmypqomFrmOx8bhMzeMyNY/L1ObjuCxEU88dwFFD4w8HvrkAdUhSkxYi+K8Gjk7yyCQHAS8CLlmAOiSpSWNv6qmqXUleC3wcuC/wvqr68rjr6NmSbKIaA4/L1DwuU/O4TG2fj8vYv9yVJC0s79yVpMYY/JLUGIN/jpK8L8mOJF/ay/yXJPli9/eZJMeOu8aFMMJxOb07Jtcn2ZLkqeOucSHMdFyGljshyd3d/S77vRFeLycl+efu9XJ9kreOu8aFMMrrpTs21yf5cpJPzmb7Bv/crQdOnWb+N4BnVNUTgP9KO19UrWf647IZOLaqjgNeBbx3DDUtBuuZ/rjs7s7kTxhc+NCK9cxwXIBPVdVx3d8fj6GmxWA90xyXJA8B3g08t6p+GThzNhs3+Oeoqq4EvjPN/M9U1T91o59lcL/Cfm+E4/KDuueKggcwxc17+6OZjkvndcBfAzv6r2hxGPG4NGeE4/LvgIuq6rZu+Vm9Zgz+8Tgb+NhCF7FYJHlekpuAv2Vw1t+8JD8PPA84b6FrWYSekuQLST6W5JcXuphF4jHAoUmuSHJNkpfPZmV/iKVnSU5mEPxNtGWPoqouBi5O8nQGzWDPWuCSFoN3AW+uqruTqXo1ada1DPqb+UGS5wAfBY5e2JIWhQOAJwGrgUOAq5J8tqq+OurK6kmSJzBow352Vd250PUsNlV1ZZJfSHJ4VbXeGdcE8KEu9A8HnpNkV1V9dEGrWmBV9b2hx3+X5N2+XoBB1zffrqofAj9MciVwLDBS8NvU05MkK4GLgJeN+i7cgiSPTpduSZ4IHAQ0/6ZYVY+sqlVVtQq4EHhN66EPkOTIodfLkxlkVvOvF2Aj8LQkByS5P/ArwI2jruwZ/xwlOR84CTg8yR3AfwEOBKiq84C3Ag8F3t29bne10NPgCMfl3wIvT3IX8GPgrKEve/dbIxyXJo1wXF4A/E6SXQxeLy/y9VLnVdWNSS4Dvgj8DHhvVU17qfC9tt/AMZQkDbGpR5IaY/BLUmMMfklqjMEvSY0x+CWpMQZ/Q7rbu39j0rTXJ3n3NOvckuTw/qubXpKndb0QXp/kkEnzpuzJMMnbk9zU9QZ6cdex1VTbXp/kG123AF9N8oGuC4Wplj0tyXXdsl9J8lsz1P2KJH8+y6c775L89u7b+ruaHjY0b8Z/46FeMq/rjuk75qmu1yb5WpJaDK+zVhj8bTmfwW8cD3tRN32xewnwjq6Hxh9PmreeqXsy3AQ8rush9avAW6bZ/puq6ljgGOA64PIMfhP6XyU5kEEvq/+mW/Z44Io5PJeRdT127rPu2u8PdKOvAB42zeJ786mqOp7B8z4tyYnzUNo/MOiy49Z52JZGZPC35UIG/2HvB5BkFYMA+HSSFye5IcmXkvzJ5BWTrBo+o07yxiR/1D2+IsmfJrkyyY0Z9Cl/UZKtSf7b0DovTfL57qz9f08VaklWd2eVN3Rn8vdL8mrghcBbk3xw8jp768mwqj5RVbu60ZF6SK2BPwW+BTx70uwHMbjp8c5u2Z9U1c1d3euTnJfkU92nhtOG1ntYksu64/E/h57rrye5Ksm1ST6S5IHd9FuSvDXJp4Ez97bc0HaWJ7mme3xsd/a8shv/xyT3T/JH3b/ZCxh0D/HBSZ+eXtdt/4YkvzjDMfoxcD0w5aei2aiq66rqln3djmbH4G9I11/Q57nn7PhFwIeBFQz6gX8mcBxwQpIzZrn5n1bV0xn0LrkR+F3gccArkjw0yWOBs4ATu77472ZwFv+vkhzM4Oz9rKp6PIOQ/Z2qei9wCYOz8nutMwuvYnY9pF4L3CsAq+o7XR23Jjk/gx/bGf4/tAp4BvCbwHnd84HBMT0LeDxwVpKjumaN/ww8q6qeCGwB3jC0rX+pqqcCfz/Dcru75D04yYOBp3XLPC3JI4AdVfWjoWUv7Oa/ZNKnp293238P8MbpDkySQxl0lHblFPOOyT0/mjL57yHTbVfjY5cN7dnd3LOxG74KOAG4oqp2AnRn1U9n0BPiqC7phjcAX66qbd22vg4cxaB30icBV2fQhcUh7Nnv/DHAN4b6NtrA4A3kXbOoYw9J/hDYBezxaWG61aaaWFWvTvJ4Bs0TbwROYdB0AnBBVf0M2No9791vHJur6p+7Wr4CPAJ4CPBLwD90x+Mg4KqhXX24G/7qDMvt9hngRAb/bv+dwZt7gE+N+Hwv6obXAM/fyzJPS/JFBv9O51bVtyYv0H0COm7EfWqBGPzt+Sjwzgw6SDukqq7d3Swwg13c+xPiwZPm/6Qb/mzo8e7xAxiE0Iaqmq6dfd77I06yBjgNWL27j5ck72fQTv3NqnrOXlY9nsGvhe2hqm4AbkjyVwx+ae0Vu2dNXrQbDh+Pu7nneGyqqhfvZf8/3P0UZlhut08xONt/BIM39Td3+790hvV2213j7vqm3EdVnZbkMQyaBy+uquuHF0hyDPe8aU12UlV9d8R61CObehpTVT9g8IXk+7jnS93PAc9IcnjX7v5iYPJveG4HlnfNNvdjEKazsRl4QZLlAEkO65oiht0ErEry6G78ZVPUMbIkpzIIwOdOau54ZdfMsUfoZ+D3GDR/XTZp3gOTnDQ06Tju/aXkmUnuk+QXgEcBN09T3meBE3c/164d/jH7sNyVwEuBrd2nju8Az2Hw5elk32fwfcWcdJ/I/geDYzt53s1DP5M4+e+7c92n5pfB36bzGfTd/SGArlnmLcDlwBeAa6tq4/AKVXUX8McM3iQuZRDSI6uqrzBoq/5E11ywiUG4Di/zL8ArgY8kuYHBp4UZe67MoCfDq4BjktyR5Oxu1p8zCLhNXRvzdNt6e5IvMLj65wTg5Kr66eRdAX+Q5OYk1wNv456zfRgE/ScZfJfw293zmVLXrPYK4PzueHyWSd8pzHK5W7qHu9vdPw18d+jnP4etZ/AdxB6Xxs7CecDTkzxyjusDkOT3Muh98uHAF5O08hvMC8reOaV5kGQ9cGn35am0qHnGL0mN8YxfkhrjGb8kNcbgl6TGGPyS1BiDX5IaY/BLUmP+PxFDpBHFVvC3AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# A natural question is: What is the uncertainty in this value? Well, we\n", "# can run this d = 12 simulation many times and plot a distribution of the\n", "# calculated prefactors. Below, we run the same simulation 400 times.\n", "d = 12\n", "n = int(1e6)\n", "maxk = 400\n", "from datetime import datetime\n", "print(datetime.now())\n", "V12List = []\n", "for k in range(maxk):\n", " Zn = 0\n", " for i in range(n):\n", " x = np.random.uniform(0, 1, d)\n", " Rsq = 0\n", " for j in range(d):\n", " Rsq = Rsq + x[j]**2\n", " R = np.sqrt(Rsq)\n", " if R <= 1:\n", " Zn += 1\n", " if k % 40 == 0: # This if statement is used to print the value of k\n", " # every 40th iteration. It helps track the progress of the loop\n", " # during execution.\n", " print(int(k/4), '%', sep = '')\n", " V12List = V12List + [2**d*Zn/n]\n", "print(datetime.now())\n", "plt.figure()\n", "nbins = 20\n", "plt.hist(V12List, nbins, color = 'c', edgecolor = 'k')\n", "plt.xlabel('Volume of 12-D Sphere with R = 1')\n", "plt.ylabel('Counts')" ] }, { "cell_type": "code", "execution_count": 15, "id": "informal-beach", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.06729470100705905" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# An estimate of the error in a Monte Carlo simulation of the volume of a\n", "# sphere in 12-D is given by the standard deviation of the V12List\n", "# distribution.\n", "import statistics as stats\n", "std = stats.stdev(V12List)\n", "std" ] }, { "cell_type": "code", "execution_count": 16, "id": "liquid-donna", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "12-D sphere volume prefactor: 1.32775936\n", "The standard deviation: 0.06729470100705905\n" ] } ], "source": [ "# However, we can do better! A better estimate of 12-D sphere prefactor is\n", "# given by the mean of V12List.\n", "prefactor12D = stats.mean(V12List)\n", "print('12-D sphere volume prefactor:', prefactor12D)\n", "print('The standard deviation:', std)" ] }, { "cell_type": "code", "execution_count": 17, "id": "filled-estonia", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The standard error: 0.0033647350503529525\n" ] } ], "source": [ "# If you've taken PHYS 232, you know that the error in the mean of a\n", "# distribution is given by the standard deviation (sigma) divided the\n", "# square root of number of trials used to produce the distribution\n", "# (N = maxk = 400 in our case). This is called the standard error.\n", "stdError = std/np.sqrt(maxk)\n", "print('The standard error:', stdError)" ] }, { "cell_type": "code", "execution_count": 18, "id": "corresponding-atmosphere", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The prefactor for the volume of a 12-D sphere is 1.32775936 ± 0.0033647350503529525\n" ] } ], "source": [ "# So we should be comfortable reporting that the prefactor used to\n", "# calculate the volume of a 12-D sphere is:\n", "print('The prefactor for the volume of a 12-D sphere is', prefactor12D,\\\n", " '\\u00B1', stdError)\n", "# Thus, we've calculated this prefactor with a prefactor with a precision\n", "# better than 1%." ] }, { "cell_type": "code", "execution_count": 19, "id": "conditional-customs", "metadata": {}, "outputs": [], "source": [ "# If you've looked at the \"Hit & Miss\" Python tutorial, you would have seen\n", "# that the error in a Hit & Miss integration decreases as 1/sqrt(N). We\n", "# first did 1e6 iterations of our calculation of the volume of a 12-D sphere.\n", "# For this calculation the error is given by the standard deviation of the\n", "# distribution that we calculated. However, in calculating the\n", "# distribution we actually ran a total of 400 million trials of our\n", "# simulation (400*1e6). Therefore, with a factor of 400 more trials, we\n", "# should expect the error in our final volume (the mean of the\n", "# distribution) to go down by a factor of sqrt(400) = 20. This is exactly\n", "# what happened when we calculated the standard error of our distribution.\n", "# Everything is beautifully consistent!" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.8" } }, "nbformat": 4, "nbformat_minor": 5 }